/**
 * @file   Interp.h
 * @author xingyifan
 * @date   2021/10/25 10:56:17
 * 
 * @brief  struct InterpCondition, class Polynomial, class NewtonInterp
 * 
 * 
 */

#ifndef __INTERP__
#define __INTERP__

#include "Config.h"
#include "Vec.h"

using namespace std;

//关于Vec的宏
#define Vec1 Vec<double, 1>
#define Vec2 Vec<double, 2>
#define VEC(Dim) Vec<double, (Dim)>

enum SplineType
{
    ppForm,
    cardinalB
};

enum BCType
{
    complete,
    notAknot,
    periodic
};

typedef struct InterpConditions
{
    vector<double> points;                   //插值点（调整为只记录端点，不重复记点）
    map<double, vector<double>> interpSites; //插值点对应函数值
} INTERPC;

template <int Order, class CoefType>
class Polynomial
{
private:
    vector<CoefType> coefficients; //系数向量

public:
    Polynomial() : coefficients(vector<CoefType>(Order + 1)) {}
    Polynomial(const vector<CoefType> &_coefficients) : Polynomial()
    {
        int size = MIN(Order + 1, _coefficients.size());
        for (int i = 0; i < size; i++)
        {
            coefficients[i] = _coefficients[i];
        }
    }
    template <int Order2>
    Polynomial(const Polynomial<Order2, CoefType> &poly) : Polynomial()
    {
        int size = MIN(Order + 1, Order2 + 1);
        for (int i = 0; i < size; i++)
        {
            coefficients[i] = poly.get_coefficients()[i];
        }
    }

//多项式加法
#define INTERP_BINARY_OP(OpName, Op)                                             \
    template <int Order2>                                                        \
    auto OpName(const Polynomial<Order2, CoefType> &poly) const                  \
    {                                                                            \
        int index = 0;                                                           \
        int this_size = order() + 1;                                             \
        int that_size = poly.order() + 1;                                        \
        vector<CoefType> coef(this_size);                                        \
        while (index < MIN(this_size, that_size))                                \
        {                                                                        \
            coef[index] = coefficients[index] Op poly.get_coefficients()[index]; \
            index++;                                                             \
        }                                                                        \
        while (index < that_size)                                                \
        {                                                                        \
            int first = 0;                                                       \
            first = first Op 1;                                                  \
            coef.push_back(poly.get_coefficients()[index] * first);              \
            index++;                                                             \
        }                                                                        \
        return Polynomial<MAX(Order, Order2), CoefType>(coef);                   \
    }

    INTERP_BINARY_OP(operator+, +)
    INTERP_BINARY_OP(operator-, -)
#undef INTERP_BINARY_OP

//多项式乘系数
#define INTERP_NUM_OP(OpName, Op)           \
    Polynomial OpName(CoefType k) const     \
    {                                       \
        vector<CoefType> coef(Order + 1);   \
        for (int i = 0; i < Order + 1; i++) \
        {                                   \
            coef[i] = coefficients[i] Op k; \
        }                                   \
        return Polynomial(coef);            \
    }

    INTERP_NUM_OP(operator*, *)
    INTERP_NUM_OP(operator/, /)
#undef INTERP_NUM_OP

    //多项式相乘
    template <int Order2>
    auto operator*(const Polynomial<Order2, CoefType> &poly)
    {
        vector<CoefType> coef;
        const int m = Order + 1;
        const int n = Order2 + 1;
        coef.resize(m + n);
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
            {
                coef[i + j] = coef[i + j] + coefficients[i] * poly.get_coefficients()[j];
            }
        }
        return Polynomial<Order + Order2, CoefType>(coef);
    }

    //不同阶多项式的赋值
    template <int Order2>
    void operator=(const Polynomial<Order2, CoefType> &poly)
    {
        int size = MIN(Order + 1, Order2 + 1);
        for (int i = 0; i < size; i++)
        {
            coefficients[i] = poly.get_coefficients()[i];
        }
    }

    bool operator==(const Polynomial &poly)
    {
        if (coefficients == poly.coefficients)
        {
            return true;
        }
        return false;
    }

    //向右平移多项式
    Polynomial<Order, CoefType> operator>>(double dx)
    {
        Polynomial<Order, CoefType> poly;
        Polynomial<Order, CoefType> origin({-1 * dx, 1});
        for (int i = 0; i < coefficients.size(); i++)
        {
            Polynomial<Order, CoefType> tmp({1});
            for (int j = 0; j < i; j++)
            {
                tmp = tmp * origin;
            }
            poly = poly + tmp * coefficients[i];
        }
        return poly;
    }

    //求函数的幂
    template <unsigned int T>
    Polynomial<Order * T, CoefType> pow()
    {
        Polynomial<Order * T, CoefType> poly(coefficients);
        Polynomial<Order * T, CoefType> origin = poly;
        for (int i = 1; i < T; i++)
        {
            Polynomial<Order * T, CoefType> tmp = poly * origin;
            poly = tmp;
        }
        return poly;
    }

    friend ostream &operator<<(ostream &out, const Polynomial &poly)
    {
        out << poly.coefficients[0];

        for (int i = 1; i < (int)poly.order() + 1; i++)
        {
            out << " + " << poly.coefficients[i] << "x^" << i;
        }
        out << endl;
        return out;
    }

    //求值
    CoefType evaluate(double x)
    {
        CoefType y = 0;
        double x_i = 1;
        for (int i = 0; i < Order + 1; i++)
        {
            y = y + coefficients[i] * x_i;
            x_i = x_i * x;
        }
        return y;
    }

    //求导
    Polynomial<Order - 1, CoefType> derivation()
    {
        vector<CoefType> _coefficients;
        for (int i = 1; i < Order + 1; i++)
        {
            _coefficients.push_back(coefficients[i] * i);
        }
        return Polynomial<Order - 1, CoefType>(_coefficients);
    }

    //不定积分
    Polynomial<Order + 1, CoefType> integration()
    {
        vector<CoefType> _coefficients(Order + 2);
        for (int i = 0; i < Order + 1; i++)
        {
            _coefficients[i + 1] = coefficients[i] / (i + 1);
        }
        return Polynomial<Order + 1, CoefType>(_coefficients);
    }

    //定积分
    CoefType integration(double left, double right)
    {
        Polynomial poly(integration());
        CoefType y = poly.caculate(left) - poly.caculate(right);
        return y;
    }

    //绘制图像
    void plot(double left, double right, char color)
    {
        //调整输出精度
        cout << setprecision(16);
        cout << "x = linspace(" << left << "," << right << ");" << endl;
        cout << "plot(x," << coefficients[0];

        for (int i = 1; i < Order + 1; i++)
        {
            cout << " + " << coefficients[i] << "*x.^" << i;
        }
        cout << ",\'" << color << "\');" << endl;
        cout << setprecision(6);
    }

    void plot2(double left, double right, char color)
    {
        cout << setprecision(16);
        cout << "t = linspace(" << left << "," << right << ");" << endl;
        cout << "plot(" << coefficients[0][0];

        for (int i = 1; i < Order + 1; i++)
        {
            cout << " + " << coefficients[i][0] << "*t.^" << i;
        }

        cout << ", " << coefficients[0][1];

        for (int i = 1; i < Order + 1; i++)
        {
            cout << " + " << coefficients[i][1] << "*t.^" << i;
        }
        cout << ",\'" << color << "\');" << endl;
        cout << setprecision(6);
    }

    const u_int order() const { return Order; }                        //获取多项式的阶
    vector<CoefType> get_coefficients() const { return coefficients; } //获取多项式系数
};

/**
 * @date   2021/11/24 11:43:50
 * 
 * @arg    const Polynomial<Order, VEC(Dim1)> &poly1, const Polynomial<Order, VEC(Dim2)> &poly2
 * @return Polynomial<Order, VEC(Dim1 + Dim2)>
 * 
 */

//将两个同阶多项式合并为高维多项式
template <int Order, int Dim1, int Dim2>
auto overlying(const Polynomial<Order, VEC(Dim1)> &poly1, const Polynomial<Order, VEC(Dim2)> &poly2)
{
    vector<VEC(Dim1 + Dim2)> coef(Order + 1);
    vector<VEC(Dim1)> coef1 = poly1.get_coefficients();
    vector<VEC(Dim2)> coef2 = poly2.get_coefficients();
    for (int i = 0; i < Order + 1; i++)
    {
        for (int j = 0; j < Dim1; j++)
        {
            coef[i][j] = coef1[i][j];
        }
        for (int j = Dim1; j < Dim1 + Dim2; j++)
        {
            coef[i][j] = coef2[i][j - Dim1];
        }
    }
    return Polynomial<Order, VEC(Dim1 + Dim2)>(coef);
}

#endif
